www.gusucode.com > 超声波测量以及形成图像 对相关信号进行模拟仿真 > 超声波测量以及形成图像 对相关信号进行模拟仿真/digital holograpy/prog/ASP.m
function varargout = ASP( U1,z,dxy1,lambda,varargin ) %ASP Diffraction calculation by Angular Spectrum Propagation method % Syntax: % [U2,dxy2] = ASP(U1,z,dxy1,'PropertyName','PropertyValue',...) % [U2,...] = ASP(U1,z,dxy1,...) % ASP(U1,z,dxy1,...) % % U1 is the wavefront of the object plane % U2 is the wavefront of the diffraction plane % U1 and U2 are all two-dimensional array % size of U1 and U2 are even % z is the distance between object plane and diffraction plane % dxy1 is the sampling distance of the object, dxy1=[dx1,dy1] % lambda is the wavelength of the laser % % if there is no output, image on diffraction plane will be displayed % else, no image will be displayed % % the origin of coordinates is at M/2+1,N/2+1 % %-------------------------------------------------------------------------- % PropertyName and PropertyValue: % % piz {off} | t % pad the input array with zeros % off - do not pad % t - paste U1 to an all zeros array whose size is t*size(U1) % % iiz {off} | t % interpolate the input array with zeros % off - do not interpolate % t - interpolate t-1 points between every two points % % pfz {off} | t % pad the frequency array with zeros % off - do not pad % t - paste the frequency array(F) to an all zeros array whose % size is t*size(F) % % FTmethod {fourier} | b,FSN,f0 % which method to be use to realize the fourier transform % fourier - traditional fourier transform using single fft % b,FSN,f0 - FAFS using 3 ffts, b,FSN,f0 are the parameters needed % % IFTmethod {invfourier} | b,SSN,xy0 % which method to be use to realize the inverse fourier transform % invfourier - traditional inverse fourier transform using single fft % b,SSN,xy0 - IFASS using 3 ffts, b,SSN,xy0 are the parameters needed %-------------------------------------------------------------------------- error(nargchk(4,14,nargin)) if nargout>2 error('Too many output arguments') end for n=1:length(varargin) if ~ischar(varargin{n}) error('Property names and values must be characters') end end % ----------------------construct the property array---------------------- parray=[0,0,0,0,0]; % initialize the property array, the property array % is used to store the property specified by the user freepv={'1','1','1','',''}; % initialize the free property value array, free property % values do not fall into those values in the property structure property=struct('name',{'piz','iiz','pfz','FTmethod','IFTmethod'}); % construct the property structure, property(1).value={'off',''}; % which stores all possible properties property(2).value={'off',''}; property(3).value={'off',''}; property(4).value={'fourier',''}; property(5).value={'invfourier',''}; for l=1:2:length(varargin) namefound=0; valuefound=0; m=0; n=0; while ~namefound && m<length(parray) m=m+1; if strcmp(varargin{l},property(m).name) namefound=1; end end while namefound && ~valuefound && n<length(property(m).value) n=n+1; if strcmp(varargin{l+1},property(m).value{n}) valuefound=1; parray(m)=n; end end %---------------------------------------------- if namefound==1 && valuefound==0 parray(m)=3; freepv{m}=varargin{l+1}; valuefound=1; end %---------------------------------------------- if namefound==0 error('wrong property name') elseif valuefound==0 error('wrong property value') end end %-------------------------------------------------------------------------- k=2*pi/lambda; if parray(1)==3 t=str2num(freepv{1}); U1=paste(zeros(t*size(U1)),U1); end if parray(2)==3 t=str2num(freepv{2}); temp=U1; U1=zeros(t*size(U1)); U1(1:t:end,1:t:end)=temp; clear temp; end %---------------------- implement fourier transform ---------------------- if parray(4)~=3 [F,dfxy(1),dfxy(2)]=fourier(U1,dxy1(1),dxy1(2)); else eval(['[F,dfxy] = FAFS(U1,dxy1,',freepv{4},');']); end clear U1 if parray(3)==3 t=str2num(freepv{3}); F=paste(zeros(t*size(F)),F); end %---------- construct the transfer function in frequency domain ---------- [M,N]=size(F); H=zeros(M,N); for m=1:M for n=1:N fx=(n-N/2-1)*dfxy(1); fy=-(m-M/2-1)*dfxy(2); H(m,n)=exp(i*2*pi*z*sqrt((1/lambda)^2-fx^2-fy^2)); end end %------------------------------------------------------------------------- F=F.*H; clear H; %------------------ implement inverse fourier transform ------------------ if parray(5)~=3 && nargout~=0 [U2,dxy2(1),dxy2(2)]=invfourier(F,dfxy(1),dfxy(2)); elseif parray(5)==3 && nargout~=0 eval(['[U2,dxy2] = IFASS(F,dfxy,',freepv{5},');']); elseif parray(5)~=3 && nargout==0 invfourier(F,dfxy(1),dfxy(2)); else eval(['IFASS(F,dfxy,',freepv{5},');']); end %-------------------------------------------------------------------------- switch nargout case 1 varargout{1}=U2; case 2 varargout{1}=U2; varargout{2}=dxy2; end